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Abstract. We present numerical simulations of the defocusing nonlinear Schrodinger (NLS) equation with 
an energy supercritical nonlinearity. These computations were motivated by recent works of Kenig-Merle 
and Kilip-Visan who considered some energy supercritical wave equations and proved that if the solution is 
a priori bounded in the critical Sobolev space (i.e. the space whose homogeneous norm is invariant under 
the scaling leaving the equation invariant), then it exists for all time and scatters. 

In this paper, we numerically investigate the boundedness of the // 2 -critical Sobolev norm for solutions 
of the NLS equation in dimension five with quintic nonlinearity. Wc find that for a class of initial conditions, 
this norm remains bounded, the solution exists for long time, and scatters. 



1. Introduction 

This paper is a numerical investigation of wellposedness and scattering properties of solutions of the 
defocusing nonlinear Schrodinger (NLS) equation, 



(1.1) hit + Ait - \u\ p - l u = 0, u:R d 



x 



in the energy supercritical setting. The notion of H s -criticality is associated with the scaling transformation 
u(x,t) — > u\(x,t) = A 2 /( p-1 'u(Aa;, X 2 t) that leaves the NLS equation invariant. We say that the problem is 
ff s -critical, if the homogeneous ff s -norm of the solution remains unchanged under the above scaling. We 
place ourselves in the energy (H 1 ) supercritical regime by assuming that the dimension and the nonlinearity 
are such that the critical Sobolev exponent s c satisfies 

(1.2) »-5-^i> 1 - 

The two conserved quantities of the NLS equation, the mass 

(1.3) M(u) = f \u{x,t)\ 2 dx 
and the energy 

(1.4) E0u)= [ (\Vu(x,t)\ 2 + \u\ p+1 )dx 

Jr* P + 1 

have indeed led to a special interest to mass-critical (s c = 0) and energy-critical (s c = 1) problems. 

Significant progress on global well-posedness and scattering has been made since the pioneering work on 
scattering by Ginibre and Velo [5], Lin and Strauss [TB], and Strauss [52]. In these early works, scattering 
was shown to hold for a range of energy subcritical configurations with finite variance. Subsequent work [9] 
proved scattering in the energy space. These results, and refinements, are collected in [23], [25], [4]. 

For energy subcritical problems with s c < 1, global well-posedness and scattering results have been 
obtained for H s data with s near 1 and away from s c . Bourgain [1 established scattering in H s , for radially 
symmetric data, for all s € (11/3, 1). This was refined by Colliander, Keel, Staffilani, Tao and Takaoka [S], 
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where the solution emerging from general initial data was shown to scatter for s € (4/5, 1). These results 
establish global well-posedness and scattering for some data which have infinite energy. 

Scattering results at critical regularity have been established starting ten years ago. A breakthrough result 
[2] on the 3d energy critical problem for radial data established a new strategy for proving scattering for 
initial data of critical regularity. With other ideas in [5], Bourgain's induction-on-energy strategy resolved 
the case for general data. The 4 dimensional defocusing energy critical case was then established by Ryckman 
and Visan [2U] with the higher dimensional case thereafter [57]. Other important advances on the energy 
critical problem appear in [10], [26]. Scattering was obtained [11] beneath the natural threshold size for the 
focusing energy critical problem. Moreover, ideas in [111 113] simplified the implementation of the energy 
critical strategy leading to advances on other model equations. Building on these developments, scattering 
in the L 2 -critical problem for large radial L 2 data was established in the defocusing case [TB] and for radial 
data under the ground state mass in the focusing case [17] . 

Scattering is expected to hold for general large data of critical regularity for ( |1.1[ ). To date, no such result 
has been established except in the energy critical case. Under the assumption of bounded critical norm, the 
defocusing H 1 / 2 case has been shown [T5] to scatter for large critical data. 

The hypothesis of a bounded critical Sobolev norm has recently been considered in the energy supercritical 
regime s c > 1. Following a recent work by Kenig and Merle [2] on the energy supercritical wav^J equation, 
Killip and Visan |16j considered some classes of energy supercritical NLS equations in dimension d > 5 and 
proved that if the solution is a priori bounded in the critical Sobolev space H Sc , it exits for all time and 
scatters. 

The purpose of our study is to investigate numerically the latter assumption on the critical Sobolev norm. 
We are also motivated by the discussion of the "theoretical possibility for computer assisted proofs of global 
well-posedness and scattering" appearing in [3] • We have performed our computations in the case d = 5 with 
quintic nonlinearity (p — 5). This is the 'simplest' cas^Jwhen the critical exponent is the smallest possible 
integer; s c — 2. In addition, we assume the initial conditions are spherically symmetric to simplify the 
computations. Note that our choice of dimension and nonlinearity is not exactly covered by Theorems 1.4 
and 1.5 of |16j . although the authors claim that their result can be extended to any power law nonlinearity 
with p odd integer, and s c < p. 

The main observation of our numerical work is that, for the various initial conditions we have considered, 
the critical norm H 2 of the solution remains bounded for all time and that the solution scatters. As time 
evolves and the solution reaches an asymptotic state, the energy concentrates into the kinetic energy and 
the potential energy tends to zero. At the same time, the H 2 norm of the solution stabilizes to some value. 
This value may be much higher than that of the initial conditions. For a more quantitative assessment of 
the solution in the frequency space, we also calculate the norm of the solution in the Besov space B 2 oa . Let 
Pj be the Fourier multiplier operator which projects the Fourier transform of a function / onto the annulus 
2 j < |£| < 2 j+1 . The Besov B| 00 (]R' i ) is equipped with the norm 

(1.5) sup2 2 ->||P,/|| i2E . 

lr The motivation for the numerical study described in this paper applies equally well to the nonlinear wave equation. As we 
had an existing code for simulating radially symmetric NLS , we chose to study NLS. In a future publication we shall examine 
the supercritical wave equation. Similar computations were performed by Strauss and Vazquez on nonlinear Klein-Gordon 
equations [24] . 

2 We observe similar phenomena in the d = 6, p = 3 case, which is also H 2 critical. 
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We observe that the H 2 density spreads out in Fourier space and the B 2 ^ norm shrinks to small values 
as time advances under (1.1 1. The emergent small Besov norm and Strichartz refinements (such as those 
appearing in the work [TS] of Planchon and references therein) gives further evidence of scattering for solutions 
of ( 1.1 1 and provides a smallness mechanism for possibly implementing the computer assisted proof described 
in 0. 

In Section 2 we describe our numerical schemes and in Section 3 we present our main results. 

2. Numerical Methods 



2.1. Time and Space Discretization. We study ( 1.1 ) with d = 5 and p = 5 and radially symmetric data. 
This configuration is convenient because the H 2 norm is bounded by the L 2 norm of Au. 

To simulate the problem, we first truncate the domain, restricting r G (0, R max ), with boundary conditions 

(2.1) u r \ r=0 =0, u(r = i? max ) = 0. 

These are interpreted as a symmetry condition at the origin, and an infinite barrier at r = i? max so that 
u = at all r > i? max - i?max must be sufficiently large to avoid boundary interaction. The domain is 
discretized as 

(2.2) = r < n < r 2 < ■ . ■ < rj = jh < . . . < r N = L. 

Letting Uj(t) — u(rj,t), the fourth order spatial discretization is 

• -Uj-2 + 16f7j_i - 30^ + mUj+i - Vj+2 
j I2h 2 
(2 ' 3) , 4 Uj-2 - Wj-i + 8U J+1 - U j+2 



rj I2h 



= \U4PU4. 



and the discretized boundary conditions are: 

(2.4) U-! = Ui ; C/_ 2 - U 2 , 

(2.5) U m+i = ; U N+2 = 0. 

We integrate in time using the classical fourth order Runge-Kutta scheme with At = 0(Ar 2 ) to ensure 
numerical stability. 

2.2. Fourier Transform. We need to compute the Fourier transform of the solution to evaluate its norm 
in the Besov space. The Fourier transform of a radial function in E rf at k = |k| is 

1 f°° d — 9 

(2.6) fi(fc) = — / u{r)J v {kr)r d/2 dr, v = 



k» J Q * — '■ ~" 2 ■ 

This formula is derived in many texts on the Fourier transform, including Stein and Weiss |2lj^[ Indeed, for 
radially symetric functions, 

The kernel is 

e -ikrl-y dSy = (2n) d ' 2 (kr)^ d - 2 y 2 J {d _ 2)/2 (kr) 1 



aso(i) 



leading to expression |2.6 



■^Stein and Weiss used a different definition of the Fourier Transform. Thus, their formula is slightly different. 
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Following Cree and Bones [7J, we approximate the integral as follows. For k > 0, we use the trapezoidal 
rule at the grid points Tj — jAr, Ar = R ma , x /N and j = 0, 1, . . . , N. At k = 0, we instead approximate the 
integral 

1 f°° 

(2.7) u(k = 0) = -, c / u(r)r v+d ' 2 dr. 

1 ' 1 ' 2T(1 + v) J Q W 

The Nyquist frequency, i^ max , is related to our discretization by the expression: 

N 

The Fourier transform is computed at the points kj = jAk, Ak = l/(2i? max ) and j = 0, 1, . . . , N. As 
discussed in their article, Cree and Bones found this method to be robust, though it is slow. 

2.3. Besov Approximation. We now approximate the Besov space norm 

(2.8) \\u\\g 2 = sup2 2j ||M|| i 2 ([23i 2 3+ i )) 

For this purpose, we first identify values of j for which these integrals can be meaningfully computed by the 
trapezoidal rule. Let 

(2.9) j min = ceil(log(4fc 1 )/log(2)) 

(2.10) j max = floor(log(fc 7V )/log(2)) 

We choose j m j„ to guarantee at least four grid points < 2- 7min . For j m i n < j < j ma x ~ L the integral 

(2.11) INIi 2([23 , 23+1)) - / H 2 dfc~<7- 
is computed by the trapezoidal rule. We also compute 

(2-12) NI^ ([ o )2 W)) = / 

and 



2J 



(I 



( 2 - 13 ) \\u\\l H[ 2^,K w] ) = i " \u\ 2 dk^q jn 

With these integrals in hand, 



(2.14) \\u\\ m «. max {2^'^} 

2.4. Error of Numerical Scheme. We have tested our scheme by varying both the domain size and the 
grid resolution. In Table [T] we show two metrics of our simulations, the value of \u\ at the origin, and the 
maximum of |u| at a fixed time. We see consistency amongst the simulations for the different parameters. 
Examining the tails of |«| in Figure [l] we get a qualitative assessment of how these parameters influence the 
simulation. So long as we have not reached the edge of the domain, which none of these simulations have, 
i? max matters little. In the better resolved simulations, |u| has propogated farther too the right. This is 
expected since greater resolution resolves higher wave numbers. As will be argued in the next section, u is 
scattering, and thus obeys the linear equation, where higher frequencies propagate with greater speed. 

The discretization scheme is neither mass nor energy conservative. A calculation of the relative error on 
these conserved quantities is another measure of accuracy. Table [2] shows that the spatial resolution for our 
various initial conditions and the relative error for the two invariants. As expected, the simulations with 
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No. Points 


■Rniax 


|u|(r = 0) 


max re r fl id 


10000 


100 


0.7126025579 


2.665689301 


20000 


100 


0.712561663 


2.665668567 


40000 


100 


0.7125588732 


2.665667313 


20000 


200 


0.7126025586 


2.665689301 


40000 


200 


0.7125616583 


2.665668567 


200000 


2000 


0.7126025579 


2.665689301 


;ence of \u\(r 


= 0) and max r |u| for Gaussian data uq = 



10e" 



at t = .02. 



Initial Condition 


No. Points 


-Rmax 


T 

± max 


max 4 |%AMass| 


max t %AEnergy 


Gaussian, u — lOe -1 " 2 


40001 


100 


0.04 


6.55478e-09 


2.67077e-08 


Gaussian, uq — 10e _r 


200001 


2000 


3.2 


1.97657e-06 


7.94582e-06 


Ring, m = 8r 2 e" r2 


32001 


100 


0.2 


2.33731e-10 


2.19067e-09 


Ring, uq = 8r 2 e _r2 


120001 


2400 


9.0 


9.61835e-07 


8.70909e-05 


Osc. Gaussian, u = 4e~ 10lr2 e~ r2 


40001 


100 


0.1 


1.35194e-08 


1.8279e-08 


Osc. Gaussian, uq = 4e~ 10lr2 e~ r2 


200001 


1000 


1.0 


2.15647e-07 


2.91404e-07 



Table 2. Error in the conserved quantities for each simulation. 



better spatial resolution, hence better temporal resolution, better conserve the invariants. These invariants 
were computed by Simpson's method on the interval [0, i? max ] with the discrete densities 

[(W,) 2 + (3t/,) 2 ] rf 

and 



-(RE/iAdiBcKl/i) - (3C/ 4 A disc .3C/,) + ((m,) 2 + (3C/,) 2 ) 6 

p + 1 



r 4 . 



P 

Adisc. is the discrete Laplacian from ( |2.3[ ). 

We also verify our time stepping and Fourier approximation by simulating a linear problem, computing 
the approximate Fourier transform, and observing that it does not change in time; see Figure [2] 

3. Numerical Observations 
In this section we present and discuss our numerical simulations for energy supercritical dcfocusing NLS, 



( 1.1 ). Throughout, our initial conditions are radially symmetric, simplifying the computations. We speculate 
that the dynamics persist for general data and for other energy supercritical configurations. 

3.1. Initial Conditions. We consider several families of initial conditions. These are: 
Gaussians: 

(3.1) uo(r) = Ae- r2 

Under the linear flow, these are well known to spread and decay in L°° . Our simulations show the 
nonlinear distortion in the shape. 
Rings: 

(3.2) u (r) = Ar 2 e~ r2 
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■ - '•',„ 


« = 100, N = 40001 


• • Rm 


„ = 2U(J, N = 40001 


X X R„ 


„ = 20U, N = 20001 




« = 2(0), N = 200001 


i 




1 

* 




X 




X 

X 

X 




X 

X 





10-" 



Figure 1. Variation in the tails of \u\ at t = .04 for Gaussian data uq = lOe 




Oscillatory Gaussians: 

(3.3) uo(r) = Ae- Qir2 e- r2 

Under the linear flow with a > 0, such an initial condition will initially focus towards the origin, 
then relax and decay. Our simulations show that the nonlinearity arrests this focusing. 



In all cases, the amplitude A is taken "large enough" so that the nonlinear effects in (1.1) will, at least 



initially, be strong and we will be outside the small data regime where scattering is known to occur |25j . We 
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present the following cases: Gaussian data with A = 10, u = 10e~ r2 ; Ring data with A — 8, u — 8r 2 e~ r2 ; 
Oscillating Gaussian data with A = 4 and a = 10, uo = 4e~ 10lr2 e~ r2 . 

3.2. Results. In all simulations, we find that after a transient period, the solution remains smooth and 
monotonically decays in amplitude. This is evidence of global well-posedness. The L°° smallness of the 
solution is a first indication of scattering, as the nonlinearity becomes a perturbation of the linear equation. 
After this transient period, the L 6 -norm also decays monotonically. Since the energy invariant is conserved, 
the potential energy is absorbed by the kinetic term, as is expected in scattering. We also study the evolution 
of the i^ 4 -norm, which is important because, by Strichartz, 

If the flow is to become asymptotically linear, we would expect L 14 to decay cx t~ 15 / 7 , the theoretical rate 
of the linear flow. Indeed, when the simulation is run for a sufficiently long time, this is observed. Finally, 
there is the critical norm, H 2 . For numerical convenience we track ||Au||^2, which controls H 2 . Finally, our 
simulations indicate that ||Au|j^2 saturates to a finite value. Thus, we have evidence that the scale invariant 
norm is globally bounded in time. This was the fundamental a assumption in [14\ 116] . 

Let us examine the profiles from our simulations. The evolution of the Gaussian data, Uq — 10e~' r is 
plotted in Figure [3] The shape is distorted, but \u\ is monotonically decreasing in time. We can also see 

9 2 

that oscillations in the shape spread. The ring data, uq = 8r e~ r , has more complex transient dynamics. 
As shown in Figure [4] there is initially a focusing of u towards the origin. This subsequently relaxes, and u 
spreads and decays in amplitude, much like the Gaussian data. The oscillatory gaussian, u — 4e~ 10lr2 e _r ' 2 , 
is similar. It initially focuses, subsequently relaxes, and appears to scatter as in Figure [5j 

For comparison, we simulate the linear Schrodinger equation with data uq = 4e _10lr e~ r and the same 
discretization as in the nonlinear problem. Figure [6] shows much stronger focusing towards the origin than 
seen in Figure [5] The nonlinearity is arresting the rush towards the origin and keeps the amplitude orders 
of magnitude smaller. 

For a more quantitative assessment of scattering, we examine the aforementioned integrated quantities. 
In Figure[jJ(a) and (b), we plot ||Aw||£2. Figure[jJ(a) shows the rapid growth of this norm to ~ 1200, orders 
of magnitude larger than the initial value. Figure [7] (b), computed for a longer time, suggests that it has 
saturated at this value. As expected, the potential energy, /|u| 6 , vanishes. We see this in in Figures [7] (c) 
and (d), where we have plotted ||u||^6 from the same two simulations. Lastly, in Figure [7] (f) the L 14 space 
norm, after a sufficient time, begins to decay as cx i~ 15 / 7 . 

The ring data and the oscillatory gaussian are, asymptotically, very similar. The same plots appear in 
Figures §and|9] We see rapid saturation of H 2 , the decay of the potential energy, and the asymptotically 
linear decay of the L 14 norm. One notable difference is that in the oscillating gaussian simulation, the initial 
focusing causes a decrease in H 2 and an increase in the potential. Furthermore, the saturated value of 
||Au||^2 is not appreciably larger than in its initial value. These simulations suggest there are at least two 
different time scales of interest. Saturation of H 2 happen very rapidly. In contrast, the expected asymptotic 
decay of L 14 sets in at a much later time. 

3.3. Fourier and Besov. Much of the recent analytical progress on NLS used careful treatment of the 
equation in the Fourier domain. We examine the Fourier transform of our simulations for hints that might 
be applied to future analysis. The transform, plotted at various times in Figures |7j(e),[8](e), and[9](e) shows 
several features. There is an initial spreading into high wave numbers, and the support is much broader than 
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Time 


\\&u\\ L i 




0.000 


20.2791 


1 


Do / oy 


0.004 


679.386 


1 




0.008 


864.094 


1 


ZDOU i 


0.012 


1040.59 


1 


1 7/1 £Q 
1 / 4Do 


0.016 


1119.92 


1 
1 


Zol4o 


0.020 


1150.59 


1 


ZD i 6 i 


0.024 


1163.46 


1 


28262 


0.028 


1169.47 


1 


28782 


0.032 


1172.59 


1 


28787 


0.036 


1174.35 


1 


28542 


0.040 


1175.41 


1 


28225 



Table 3. Comparison of norms for the Gaussian data with i? max = 100 and N = 10000. 
Both Besov and Sobev saturate very rapidly. 



Time 


\\Au\\ L 2 




0.000 


17.6789 


1 


66075 


0.002 


43.9913 


1 


55944 


0.004 


63.3487 


1 


59434 


0.006 


74.0549 


1 


55922 


0.008 


77.8274 


1 


55706 


0.010 


79.2784 


1 


56396 


0.012 


80.0381 


1 


56898 


0.014 


80.5317 


1 


57103 


0.016 


80.8207 


1 


57138 


0.018 


80.9811 


1 


57106 


0.020 


81.0696 


1 


57057 



Table 4. Comparison of norms for the ring data with i? max = 100 and N = 32000. Both 
Besov and Sobev saturate very rapidly. 



the initial condition. This relaxes, and the limiting state has a smaller support than during the transient 
period, but still in excess of the initial condition. The asymptotic constancy of the transform is further 
evidence that u evolves linearly and scatters. 



We can also interpret the Fourier data through the Besov norm (1.5 1. Using the method described in 
we approximate the Besov norm B 2 ^ at several times in each simulation. This data, appearing 



2.3 



Section 

as x's in Figures [7] (a), [8] (a), and [9] (a), shows several things. First, the Besov norm flf )0o is orders of 
magnitudes smaller than the scale invariant norm H 2 . Like the Sobolev norm, there is some transient 
variation followed by saturation to some asymptotic value. In the case of the oscillating gaussian data, 
Figure [9] (a), the dynamics of the two norms appear to be phase locked. Another feature is that while the 
Sobolev norm can increase (substantially) as it saturates, the Besov norm always decays. More detailed 
information for each of the three simulations is given in Tables [3j |4j and [5j 
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Time 


||A U || i2 


\\ U \\ R2 


0.00 


Q 1 O 


Oil 
Li ( 


U 


fifi cool 


0.01 


Doy 


79 Q 


u 


OUOoZO 


0.02 


7GQ 


yo 


u 


004Z4o 


0.03 


oZO 


OOZ 


u 


004oOZ 


0.04 


C07 


A A O 

44 y 


n 
U 


004Z 


0.05 


G97 
oZ ( 


/ICQ 
4oo 


n 

U 


DD4ZDZ 


0.06 


R97 

OZi i 


L ±o\) 


n 
u 




0.07 


827 


487 





664269 


0.08 


827 


487 





664269 


0.09 


827 


487 





664269 


0.1 


827 


487 





664269 



Table 5. Comparison of norms for the oscillating gaussian data with R max = 100 and 
N = 40001. Both Besov and Sobev saturate very rapidly. 
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NUMERICAL SIMULATIONS OF THE ENERGY- SUPERCRITICAL NONLINEAR SCHRODINGER EQUATION 
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Figure 3. Evolution of u = lOe r . Computed on the domain [0, 100] with 40000 + 1 points. 




Figure 4. Evolution of u = 8r 2 e r . Computed on the domain [0, 100] with 32000 + 1 points. 




Figure 5. Evolution of u = 4e 10ir "e r . Computed on the domain [0, 100] with 40000 + 1 points. 
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Figure 6. Evolution of uq = 4e 10lr e r under the linear flow. Computed on the domain 
[0, 100] with 40000 + 1 points. 




Figure 7. Metrics for u = 10e r . Figures (b), (d), (f) are computed on the domain 
[0,2000], with 200000+1 points. 
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Figure 8. Evolution of uq — 8T 2 e~ r . Figures (a), (c), and (e) are computed on the domain 
[0, 100] with 32000+1 points. Figures (b), (d), and (f) are computed on the domain [0, 2400], 
with 120000+1 points. 



NUMERICAL SIMULATIONS OF THE ENERGY- SUPERCRITICAL NONLINEAR SCHRODINGER EQUATION 
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Figure 9. Evolution of uq = 4e 10lr e r . Figures (a), (c), and (e) are computed on the 
domain [0, 100] with 40000+ 1 points. Figures (b), (d), and (f) are computed on the domain 
[0, 1000], with 200000+1 points. 



